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Abstract 

Magnetic Resonance Imaging has become nowadays an indispensable 
tool with applications ranging from medicine to material science. How- 
ever, so far the physical limits of the maximum achievable experimental 
contrast were unknown. We introduce an approach based on principles 
of optimal control theory to explore these physical limits, providing a 
benchmark for numerically optimized robust pulse sequences which can 
take into account experimental imperfections. This approach is demon- 
strated experimentally using a model system of two spatially separated 
liquids corresponding to blood in its oxygenated and deoxygenated forms. 

Since its discovery in the forties, Nuclear Magnetic Resonance (NMR) has 
become a powerful tool [TJ [5] to study the state of matter in a variety of domains 
extending from biology and chemistry [3] to solid-state physics and quantum 
computing [5J [5] . The power of NMR techniques is maybe best illustrated by 
medical imaging |6J, where it is possible e.g. to produce a three-dimensional pic- 
ture of the human brain. NMR spectroscopy and Magnetic Resonance Imaging 
(MRI) involve the manipulation of nuclear spins via their interaction with mag- 
netic fields. All experiments in liquid phase can be described in a first approach 
as follows. A sample is held in a strong and uniform longitudinal magnetic 
field denoted Bq. The magnetization of the sample is then manipulated by a 
particular sequence of transverse radio-frequency magnetic pulses B\ in order 
to prepare the system in a particular state. The analysis of the radio-frequency 
signal that is subsequently emitted by the nuclear spins leads to information 
about the structure of the molecule and its spatial position. One deduces from 
this simple description that the crucial point of this process is the initial prepa- 
ration of the sample, i.e. to design a corresponding pulse sequence to reach 
this particular state with maximum efficiency. The maximum achievable effi- 
ciency can be determined for the transfer between well defined initial and target 
states [7] if relaxation effects can be neglected. In imaging applications, where 
relaxation forms the basis for contrast, a very large number of different strate- 
gies have been proposed and implemented so far with the rapid improvement 
of NMR and MRI technology [6]. However, there was no general approach 
to provide the maximum possible performance and the majority of these pulse 
sequences have been built on the basis of intuitive and qualitative reasonings 
or on inversion methods such as the Shinnar-Le Roux algorithm [S]. Note that 
this latter can be applied only in the case where there is no relaxation effect 
and radio-frequency inhomogeneity. 

A completely different point of view emerges if this problem is approached 
from an optimal control perspective. Optimal control theory was created in its 
modern version at the end of the 1950s with the Pontryagin Maximum Principle 
(PMP) [§J[Tni[TT]. Developed originally for problems in space mechanics, opti- 
mal control has become a key tool in a large spectrum of applications including 
engineering, biology and economics. Solving an optimal control problem leads 
to the determination of a particular trajectory, that is a solution of an associated 
Hamiltonian system constructed from the PMP and satisfying given boundary 
conditions. This approach has found remarkable applications in quantum com- 
puting and NMR spectroscopy, but its application to MRI has been limited to 
the numerical design of slice-selective 90° and 180° pulses [12J. 

Despite the efficiency of MRI techniques currently used in clinics, some as- 
pects still pose fundamental problems of both theoretical and practical interest. 
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The enhancement of contrast remains one of the crucial questions for improving 
image quality and the corresponding medical diagnosis. The use of particular 
pulse sequences to generate image contrast based on relaxation rates is not new 
in MRI, since this question was raised at the beginning of the development of 
MM in the 1970s. Different strategies have been proposed, such as the Inversion 
Recovery sequence [131 [2] for T\ contrast and pulses for ultra short echo time 
experiments for T 2 contrast [T5] (See Eq. fT]) for the definition of T\ and T2 
parameters). Here, we go beyond such intuitive methods by using the powerful 
machinery of optimal control, which provides in this case not just an improved 
performance but an estimate of the global optimum, i.e. the best possible con- 
trast within the experimental constraints (see the supplementary material for 
mathematical details) . This optimized contrast is demonstrated in a laboratory 
benchmark experiment. 

In its simplest form, the contrast problem can be stated by assuming that the 
signal is composed of two different contributions. We consider as a benchmark 
example the case of (a) oxygenated vs. (b) deoxygenated blood, where the 
spins that are probed are the ones of the hydrogen atoms of water (H2O). This 
is e.g. an important issue in functional studies of the human brain. These spins 
have different relaxation rates due to the interaction with other molecules such 
as hemoglobin in its oxygenated or non-oxygenated form, leading thus to two 
different signatures of the relaxation dynamics of the magnetization, which is 
governed by the Bloch equations: 

-ojAP y + oj y AP z -M x /T* 

uM 1 x -uj x M 1 z - M l ylT l 2 (1) 
-u v M* + u x M l y - {Ml - M )/Tf , 

where M % = (M*,M*,M*) is the magnetization vector considered with % = 
(a, 6), Mq is the equilibrium magnetization, T[ and T% the longitudinal and 
transverse relaxation rates, u> the resonance offset and oj x and oj y the components 
in a rotating frame of the transverse magnetic field along the x- and y- directions. 
While different definitions of contrast exist in the literature [5] , here we consider 
a particular case that we call the saturation contrast, where the objective of 
the control problem is to find the pulse sequence which completely suppresses 
the contribution of one of the two magnetization vectors while maximizing the 
modulus of the other. If such a pulse module can be found, it can be used in 
combination with a large number of possible host sequences for imaging and 
spectroscopy (see Chapter 14 of Ref. [6] for details). 
Results 

In this section, we first analyze the ideal situation of a homogeneous ensemble 
of spin 1/2 particles irradiated on resonance, which is described by Eq. (TTJ) with 
uj = 0. We introduce the normalized vectors V 1 — M l /Mq with coordinates 
(X 1 , Y l , Z l ) to eliminate the equilibrium magnetization Mq. Based on numerical 
computations (see also [TOJ [T7J [T5] for analytical details), we restrict, without 
loss of generality, the dynamics to a meridian plane by assuming u> y — 0. Using 
advanced techniques of geometric optimal control theory [TUl [T5] , it is possible 
to find the desired control field ui x (t) which provides the optimum contrast by 
a direct solution of the PMP. The details of the theoretical approach are given 
in Section 1 of the supplementary material. 



dt 

dt 
dMl 
dt 
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As a test case, we chose the typical relaxation parameters for (a) oxygenated 
and (b) deoxygenated blood with identical T\ values (T"' 6 = 1.35 s) but dif- 
ferent T2 values (T-f = 200 ms, T| = 50 ms). Note that in this situation, the 
conventional Inversion Recovery experiment, which relies on T\ differences can- 
not provide useful contrast. However, applying the optimal control approach for 
these parameters the maximum of the modulus of V a (under the condition that 
\V h \ = 0) is found to be \V a \ = 0.4663, representing the maximum achievable 
saturation contrast in this system. Note that a similar computation can be done 
to saturate the spin (a) and maximize \V b \, with a final result of \V b \ = 0.4731. 
We underline that this saturation contrast is optimal in the sense that it is the 
physical upper limit that can be reached within the experimental constraints 
given here by the values of the relaxation rates. The shape of the optimal pulse 
is shown in Fig. [IJi. For our demonstration experiments, we did not use actual 
blood samples but prepared two different solutions with similar physical charac- 
teristics (see the Methods section). As the experimental T 2 * values [2] of the test 
sample were only about 2/3 of the T2 values assumed for the optimized pulse, 
the pulse duration and amplitude were scaled by 2/3 and 3/2, respectively, re- 
sulting in a scaled pulse duration of 0.214 s and a maximum pulse amplitude of 
the order of 10 Hz (see Fig. QJi). Within the experimental accuracy, we have 
checked that, in this case, this scaling procedure provides the same contrast as 
the optimal solution. The optimal control held was implemented experimen- 
tally as a shaped pulse on a standard Bruker Avance III 600 MHz spectrometer 
and the experimental trajectories were measured using two different samples 
in different test tubes in order to approach ideal experimental conditions with 
negligible magnetic held inhomogeneities. Under these conditions, the dynamics 
is described with high accuracy by the Bloch equations ([T]). The simulated and 
experimental trajectories of the two magnetization vectors V a (t) and V h (t) are 
shown in Fig. [TJd and a good match is found between theory and experiment. 

So far, we assumed that there is no experimental imperfection due to mag- 
netic held inhomogeneities, which are however not negligible in realistic imaging 
experiments and have therefore to be taken into account. Nevertheless, it is im- 
portant to point out that the analytical PMP-based optimal solution in the 
absence of experimental imperfections gives an estimate of a previously unavail- 
able physical upper limit of the maximum achievable saturation contrast, the 
inhomogeneities having a detrimental effect on the final result. This bound, 
which can be determined for any set of relaxation parameters, is thus relevant 
for any imaging contrast problem. Hence, it gives a fundamental benchmark to 
assess the performance of standard methods and numerically optimized pulse 
sequences in this domain. 

To perform a realistic imaging experiment, we designed a test sample consist- 
ing of a small test tube in a larger tube with an outer diameter of 8 mm, forming 
two compartments filled with solutions (a) and (b) corresponding to the relax- 
ation rates of deoxygenated and oxygenated blood, see Fig. [5] for a schematic 
representation. The experiments were performed using the same spectrometer 
as described above, equipped with a micro-imaging unit. The details are given 
in Section 2 of the supplemental material. Figure [3] shows the experimentally 
measured spatial Bq and B\ distributions in the central slice of the sample. The 
variation of Bq corresponds to resonance frequency shifts ui between and -30 
Hz, while the experimental scaling of the B\ field (which is proportional to the 
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Figure 1: (Color online) Optimal pulse sequence and trajectories for 
negligible Bq and B\ inhomogeneities. (a), The control amplitude u x (t) = 
tu x (t) /(2n) is shown for the optimal pulse sequence to maximize \V a \ under 
the condition that \V b \ — 0. (b), Corresponding simulated (solid curves) and 
experimental (open diamonds) trajectories of V a (t) (red) and V b (t) (blue) in 
the (Y, Z) plane for a homogeneous ensemble of spin particles. 




Figure 2: (Color online) Geometry of the test sample used for the imag- 
ing experiments. In the micro imaging experiments, the sample consists of 
two test tubes with outer diameters of 5 mm and 8 mm. The outer and inner 
volumes were respectively filled with the two solutions corresponding to (a) oxy- 
genated and (b) deoxygenated blood. The two slices represent the experimental 
results after the saturation of the samples (a) (top) and (b) (bottom), (see also 
the results of Fig. @| 
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Figure 3: (Color online) Experimental spatial Bq and B\ distributions. 

Spatial distributions of the Bq (a) and B\ (b) amplitudes in the central slice of 
the sample. The Bq variation is represented by the corresponding 1 H frequency 
offsets in units of Hz as a function of the spatial position, while the B\ variation 
is described by a dimensionless scaling factor. 



control amplitude ) is ±20%. 

In general, the frequency offsets created by the Bo inhomogeneities are neg- 
ligible if they are dominated by the amplitude of the control field (in units of 
Hz). However, the control amplitude of the analytically optimal pulse sequence 
shown in Fig. is less than 10 Hz, which is smaller than the resonance offset 
variation due to Bq inhomogeneities (see Fig. [3^) and therefore, the optimal 
sequence derived analytically for an ideal case is not expected to work in the 
experimental micro imaging setting. In addition, the experimental variation of 
B\ scaling was not considered in the analytical solution and is also expected to 
have detrimental effects on the pulse performance. 

In order to take into account the experimentally measured Bq and B\ distri- 
butions, numerically optimized pulse sequences were computed with the GRAPE 
algorithm [20J, which is a standard iterative algorithm to solve the optimiza- 
tion equations. The pulses are designed to work for an ensemble of spins ap- 
proximately within the range of the Bq and B\ inhomogeneities experimentally 
measured. In the numerical optimizations, we fixed the pulse duration to the 
duration of the corresponding analytical pulse. Compared to the fundamental 
contrast benchmark (|V a | = 0.47 and \V h \ = 0) provided by the PMP- based 
analytical solution, the minimum value (worst case) of \V a \ for the considered 
range of Bq and B\ inhomogeneities is 0.37 (that is 79% of the physically max- 
imum saturation contrast achievable) while the maximum value (worst case) of 
the incompletely suppressed \y h \ is 0.054. Conversely, when the goal is to satu- 
rate spin (a) and to maximize the magnetization vector of spin (b), the optimal 
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Figure 4: (Color online) Experimental implementation of the robust op- 
timal pulse. Optimization of the contrast when the deoxygenated (b) or the 
oxygenated (c) blood is saturated. Figure (a) displays the reference image after 
the application of a 90 degree pulse. 
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pulse sequence yields \V a \ = 0.047 and \V b \ = 0.33. Figure 0] shows resulting 
experimental images which are in good agreement with simulated data (see the 
supplementary material for the details of the computation). 
Discussion 

We demonstrated to which extend saturation contrast based on different relax- 
ation times T\ and T 2 can be maximized in magnetic resonance imaging within 
given experimental constraints. Starting from the analytic optimal solution of 
the homogeneous case, we have then designed a particular pulse sequence using 
numerical tools of optimal control theory to approach in a realistic experiment 
this physical limit. We emphasize that one of the main advantages of this con- 
trast enhancement is its general character since the optimal control fields can be 
computed with standard routines published in the literature and implemented 
on a standard NMR spectrometer without requiring specific materials and pro- 
cess techniques. 

The efficiency of this approach was shown in a laboratory experiment using 
a model system for the relaxation parameters of deoxygenated and oxygenated 
blood. The presented method fully exploits the differences of both T\ and 
T2 to create the maximum possible saturation contrast as opposed to conven- 
tional approaches based on T\ or T 2 differences. The combined analytical and 
numerical optimal control approach is not limited to the definition of satura- 
tion contrast (motivated by typical magnitude mode imaging experiments) used 
here for demonstration, but can also be applied to more general definitions of 
relaxation-based contrast, e.g. for phase-sensitive images and for a wide variety 
of possible host imaging sequences that can be applied after the contrast pulse 
module [B]. Furthermore, the flexibility of the optimal control approach makes 
it possible to include experimental constraints such as bounds on the control 
amplitude or pulse energy or non- linear effects such as radiation damping [21] . 

We expect that the presented principles will find practical applications in 
MRI and in particular in medical imaging, where increased contrast and sensi- 
tivity could not only help in the diagnosis but could also reduce the required 
concentration of commonly used contrast agents, which could be beneficial to 
the patient. 
Methods 

Experimental sample. The relaxation properties of oxygenated blood were 
approached by solution (a) consisting of 90% D2O, 10% glycerol and doped 
with CUSO4 with relaxation times Tf = 1.8 s, T% = 260 ms (as determined 
from CPMG-experiment) and T 2 *'° = 100 ms (as determined from the experi- 
mental line width). Deoxygenated blood was modeled by solution (b) consisting 
of 70% D 2 0, 30% glycerol and doped with CuS0 4 with T\ = 1.4 s, T| = 60 ms 
and T2 ,b =30 ms. 

Measurement method of the Bo and B\ field maps. Mapping of local Bq 
offsets was accomplished by evaluating the signal phase evolution between two 
echoes acquired in a dual-echo gradient pulse sequence f6J. The echo times of 
the 3D image acquisition were TE\ =1.5 ms and TE2 = 11.5 ms. Figure 3a of 
the main text shows the Bq field map in the central axial slice. The amplitude 
of the B\ excitation field was measured by using the cosine-like dependence of 
the remaining signal after a saturation pulse and fitting to a curve measured 
with multiple saturation flip angles |22) . The applied saturation flip angles were 
10°, 20°, . . . , 300°. Robust B\ mapping was achieved by fitting a signal model to 



the acquired slice-selective gradient echo signal in a linear least-squares sense. 
Application of the optimized pulse sequence. For the imaging experi- 
ment, the H2O content was increased in the sample in order to have a sufficient 
signal-to-noise ratio. The outer and inner volumes of the sample were filled 
with the following solutions: "oxy sample II" (80% D 2 0, 10% H 2 0, 10% glyc- 
erol doped with CUSO4 with relaxation times of Ti=2.6 s and T 2 *=100 ms) and 
"deoxy sample II" (60% D 2 0, 10% H 2 0, 30% glycerol doped with CuS0 4 with 
relaxation times of T\=1A s and T 2 =30 ms). The optimal pulses were imple- 
mented into a gradient echo pulse sequence |23| without slice selection. The 
experiments were performed with a 600MHz Bruker Avance III spectrometer 
equipped with a micro-imaging unit. The field of view (FOV) was 15mm xl5 
mm, and in the third dimension it was limited to approximately 10 mm by coil 
sensitivity. The repetition time (TR) is 7 s and the echo time (TE) is 2 ms with 
the matrix size 128x128. Figures SI and S2 of the supplementary material show 
simulated (d) and experimental (e) images using the optimized pulse to satu- 
rate the "deoxy sample II" in the inner cylinder (Fig. SI) and for the optimized 
pulse to saturate the "oxy sample II" in the outer cylinder (Fig. S2). 
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